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OBJECTIVES 

The  objectives  of  this  work  are  to  develop  and  Implement  numerical  tec hr 
nlques  that  will  enable  us  to  perform  simulations  of  spatially  evolving 
phenomena.  In  particular,  we  are  Interested  In  studying  the  phenomenon  of 
turbulent  diffusion  flame  lift-off. 

The  flame  lift-off  phenomenon  occurs  when  the  speed  of  the  fuel  jet 
exceeds  a  certain  value.  The  flame  detaches  from  the  exit  and  is  stabilized 
at  a  certain  distance  downstream.  Different  mechanisms  have  been  proposed  to 
explain  the  lift-off  phenomenon.  Generally,  it  is  believed  that  the  phenomen¬ 
on  is  similar  to  that  In  premixed  gas  combustion  (Brzustowski,  1980;  Gunther 
et  al.,  1981).  At  the  jet  exit,  the  velocity  of  the  fluid  in  the  mixing  layer 
is  higher  than  the  flame  propagation  speed.  The  flame  cannot  be  stabilized  at 
that  location.  As  the  mixing  proceeds,  the  velocity  of  the  fluid  in  the 
mixing  layer,  where  the  species  ratio  roughly  becomes  stoichiometric, 
decreases  to  the  flame  speed.  The  flame  Is  stabilized  at  that  position.  This 
model  assumes  that  the  mixing  of  the  fuel  and  the  oxidizer  as  predicted  by  the 
nonreacting  turbulence  model  reaches  the  molecular  level  everywhere  in  the 
mixing  layer.  Direct  numerical  simulations  of  a  reacting  mixing  layer  (Riley 
et  al.,  1986)  indicated  that  this  may  not  be  the  case. 

Another  model  was  suggested  by  Peters  and  Williams  (1983).  They  proposed 
that,  in  a  turbulent  mixing  layer,  the  turbulent  eddies  produce  highly 
stretched  and  contorted  sheets  across  which  molecular  diffusion  of  species 
occurs.  Combustion  appears  as  laminar  flamelets  on  these  sheets.  Based  on 
the  theory  that  a  laminar  flame  will  extinguish  itself  if  the  strain  rate  is 
sufficiently  high  that  the  local  reduced  Damkohler  number  is  lower  than  a 
critical  value  (Linan,  1974),  Peters  and  Williams  predicted  that  the  strain 
rate  near  the  exit  of  the  jet  is  too  high  for  the  flame  to  exist.  As  the 
mixing  layer  grows,  the  strain  rate  is  reduced  to  a  sufficiently  low  level 
that  a  flame  can  be  sustained.  They  proceeded  to  estimate  the  strain  rate  in 
the  mixing  layer.  Through  a  statistical  theory  that  predicts  disruption  of  a 
continuous  sheet  due  to  "boles"  on  the  sheet,  they  were  able  to  obtain 
quantitative  results  for  lift-off  distance.  Their  theoretical  predictions  of 
lift-off  heights  are  of  the  right  "order  of  magnitude"  for  limited  methane 
flame  data. 
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The  results  of  Peters  and  Williams ,  based  on  simple  turbulence  models,  are 
encouraging.  However,  a  detailed  numerical  simulation  of  the  interaction 
between  fluid  dynamical  and  combustion  processes  that  occur  in  diffusion 
flames  would  be  Invaluable  to  further  our  understanding  of  the  physics. 

Direct  numerical  simulations  consist  of  accurately  solving  appropriate 
convection-diffusion-reaction  transport  equations  by  means  of  very  accurate 
numerical  algorithms  so  that  no  turbulence  modeling  is  required.  The  direct 
numerical  simulation  technique  has  recently  been  successfully  applied  to 
chemically  reacting  flows.  Riley  et  al.  (1986)  considered  the  three- 
dimensional  temporally  growing  mixing  layer  under  the  simplest  possible 
assumption  of  a  constant-rate  chemical  reaction  with  no  heat  release.  The 
main  contribution  of  this  work  is  the  understanding  of  the  effects  of 
three-dimensional  mixing  and  diffusion  of  the  species  on  the  chemical 
reaction.  McMurtry  et  al.  (1985)  considered  the  effects  of  chemical  heat 
release  and  the  resulting  density  variation  on  the  fluid  motion  for  a 
two-dimensional  mixing  layer.  The  fluid  dynamics  and  the  chemical  reaction 
are  truly  coupled  in  this  work,  and  the  Interplay  between  the  two  are 
discussed.  However,  the  assumption  of  a  constant  chemical  reaction  is  still 
employed.  We  intend  to  extend  the  simulations  to  study  the  local  quenching 
and  the  lift-off  of  a  diffusion  flame. 

The  existing  methods  must  be  modified  to  Include  (1)  a  temperature- 
dependent  chemical  reaction  rate  and  (2)  a  spatially  developing  flow.  The 
objective  of  our  first  year's  work  is  to  examine  these  features  in  a  two- 
dimensional  mixing  layer  simulation.  The  results  of  the  present  work  can  then 
be  extended  to  three-dimensional  simulations  in  subsequent  years. 

STATUS  OP  THE  RESEARCH 

Our  efforts  in  the  past  year  were  concentrated  in  two  main  areas:  (1)  the 
incorporation  of  a  temperature-dependent  reaction  rate  into  a  previously 
developed  computer  code  and  the  eduction  of  relevant  physics  related  to  flame 
extinction  from  the  simulation  and  (2)  the  development  of  new  numerical 
algorithms  for  computing  a  spatially  developing  reacting  shear  layer,  which 
will  enable  us  to  address  the  flame  stability  and  lift-off  phenomena  in  the 
future.  In  the  first  area,  temporally  evolving  mixing  layers  were  considered, 
while  in  the  second  area,  spatially  developing  plane  jets  and  mixing  layers 
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were  considered.  In  both  areas,  relevant  experimental  measurements  were 
reviewed  and  qualitative  comparisons  with  data  were  made  whenever  possible. 

In  this  section,  a  brief  review  of  our  results  and  the  major  conclusions  drawn 
from  the  research  to  date  in  both  areas  is  given.  A  work  plan  for  the  coming 
year  is  also  presented.  Appendices  A  and  B  include  detailed  descriptions  of 
the  research. 

Summary  of  Simulations  of  Temporally  Growing  Flows 

The  mechanism  of  flame  extinction  was  studied  through  a  two-dimensional 
simulation  of  a  temporally  growing  mixing  layer  (Appendix  A).  A  binary, 
single-step,  Arrhenius-type  chemical  reaction  with  finite  rate  chemical  kinetics 
between  two  unpremixed  species  was  considered.  A  previously  developed  two- 
dimensional  numerical  code  was  modified  to  include  the  temperature-dependent 
reaction  rate  constant.  The  simulations  were  performed  with  a  fairly  large 
Zeldovlch  number  (Ze  ■  20  based  on  the  activation  energy  and  ambient  tempera¬ 
ture),  a  moderate  Reynolds  number  (Re  “  200  based  on  the  shear  layer  thickness 
and  the  velocity  difference  across  the  layer),  and  a  moderate  Damkohler  number 
(Da  »  10  based  on  the  pre-exponential  factor,  the  velocity  difference  across 
the  layer  and  the  shear  layer  thickness).  The  nonequilibrium  phenomenon 
leading  to  the  local  quenching  of  the  diffusion  flame  was  simulated  on  a 
64  X  64  grid  computational  domain.  The  results  indicate  that  the  large 
coherent  structures  have  a  major  influence  on  the  local  extinction  of  the 
diffusion  flame.  It  was  found  that  the  roll  up  of  the  unstable  shear  layer 
creates  regions  with  high  dissipation  rates  at  the  braidB  of  the  large-scale 
structure,  where  local  flame  extinction  occurs.  The  results  indicate  that  the 
local  mixture  temperature  drops  to  a  value  close  to  ambient  at  the  braids  and 
the  reaction  rate  reduces  to  zero,  even  though  the  reactants  are  well  mixed 
there.  This  is  consistent  with  the  theory  of  Peters  and  Williams  (1983). 

Within  the  limitations  of  our  numerical  resolution,  comparisons  of  the 
critical  scalar  dissipation  rate  with  that  obtained  by  large  activation  energy 
asymptotics  showed  reasonable  agreement. 

Summary  of  Simulations  of  Spatially  Developing  Flows 

In  this  task,  we  are  mainly  interested  in  the  development  of  numerical 
codes  that  are  capable  of  simulating  Bpatially  developing  flows.  In 
particular,  the  implementation  of  the  inflow/outflow  boundary  conditions  and 
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their  effect  on  the  interior  of  the  computational  domain  is  addressed. 
Presently,  we  have  finished  constructing  a  numerical  code  that  is  capable  of 
simulating  two-dimensional  mixing  layers  and  plane  jets.  For  the  purpose  of 
demonstration,  this  code  was  used  for  the  simulation  of  the  asymmetric  mixing 
in  parallel  two-stream  mixing  layers  and  was  used  to  examine  the  role  of  large 
coherent  structures  in  the  enhancement  of  chemical  reactions  that  occur  in 
such  flows  (Appendix  B). 

The  flow  in  the  cross-stream  direction  may  be  assumed  periodic.  There¬ 
fore,  pseudospectral  discretization  was  used  in  that  direction.  In  the 
8treamwlse  direction,  however,  the  flow  is  not  periodic.  Appropriate 
Inflow/outflow  boundary  conditions  must  be  applied.  At  the  inflow,  a  small 
perturbation,  in  the  form  of  the  most  unstable  mode  and  its  subharmonics,  was 
added  to  a  mean  hyperbolic  tangent  velocity  profile.  This  pertubation 
initiated  vorticlty  rollup  and  pairing  of  the  vortices  in  the  computational 
domain.  The  downstream  outflow  boundary  is  a  computational  antifact  and  is 
not  a  natural  flow  boundary.  Therefore,  the  boundary  conditions  are  difficult 
to  formulate  in  a  rigorous  sense.  So  far,  we  have  employed  a  zero-gradient 
condition  on  the  velocity  field  and  the  species  field  and  have  not  tested 
other  possibilities.  However,  computations  were  performed  with  different 
sizes  of  the  computational  domain  in  the  streamwlse  direction,  and  the  results 
indicated  that  this  artificial  boundary  condition  at  the  outflow  does  not  seem 
to  Influence  significantly  the  flow  field  in  the  interior  of  the  computational 
domain.  To  construct  the  first  working  code,  we  employed  an  overall  second- 
order  finite-difference  scheme.  The  convective  term  is  discretized  by  a 
quadratic  upwind  difference  scheme  in  which  the  leading  dispersive  truncation 
error  of  the  second-order  central  scheme  Is  removed.  This  results  in  an 
improvement  of  the  capability  in  resolving  sharp  gradients  in  a  moderate 
Reynolds  number  flow  over  the  second-order  central  scheme.  It  also  improves 
the  stability  property  over  the  latter.  Second-order  central  differencing 
requires  a  maximum  cell  Reynolds  (Peclet)  number  of  2,  which  requires  an 
excessive  and  nonpractical  number  of  grid  points  in  the  streamwlse  direction. 
With  the  upwind ing  scheme,  we  were  able  to  simulate  a  two-dimensional  flow 
with  Reynolds  and  Peclet  numbers  of  about  50  in  a  domain  with  128  X  64  grid 
points  successfully  and  with  no  numerical  oscillations.  Employing  a  spectral 
element  method  (presently  under  construction)  in  the  streamwlse  direction 
should  allow  us  to  simulate  flows  with  higher  Reynolds  numbers. 
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The  resulting  code  was  used  to  study  the  effects  of  large-scale  vortex 
dynamics  on  a  spatially  developing,  reacting  mixing  layer.  A  fast  chemistry 
I  model,  based  on  the  solution  of  a  conserved  Shvab-Zeldovich  scalar  variable, 

was  used  to  describe  the  effects  of  the  chemical  reactions.  Note  that  making 
such  an  assumption  about  the  chemical  reactions  will  not  allow  us  to  investi¬ 
gate  the  the  nonequilibrium  effects  leading  to  the  local  quenching  of  the 
§  diffusion  flame.  However,  within  this  assumption,  important  information  on 

the  convolution  of  the  flame  can  still  be  obtained.  The  results  of  the 
numerical  simulations  indicate  that  the  large  coherent  structures  have  a  major 
Influence  in  stretching  the  stream  surface  area  on  which  chemical  reactions 
0  occur,  thereby  increasing  the  total  amount  of  product  formed  in  the  reaction 

zone.  Furthermore,  the  results  exhibit  the  experimentally  observed  asymmetric 
properties  of  the  entrainment  mechanism  in  the  core  of  the  layer  (Appendix  B). 

Current  Activities 

We  are  currently  incorporating  the  temperature-dependent  reaction  rate 
into  the  code  for  spatially  developing  flow.  The  resulting  code  should  enable 
us  to  simulate  the  local  quenching  phenomena  in  such  a  flow.  We  are  also 
developing  a  three-dimensional  code  with  the  option  of  using  either  a  finite- 
difference  or  a  spectral  element  method  in  the  streamwlse  direction.  The 
resulting  code  can  be  use  to  study  three-dimensional  turbulent  fields  in  a 
spatially  evolving  mixing  layer  and  the  structure  of  "holes"  on  the  flame 
sheet.  Variable  density  effects  caused  by  the  heat  release  will  then  be  added 
to  the  code  at  a  later  stage. 

PUBLICATIONS 

The  following  written  manuscripts  have  resulted  from  our  efforts  in  the 
first  year: 

1.  Givi,  P. ,  Jou,  W.-H.,  and  Metcalfe,  R.  W. ,  "Flame  Extinction  in  a 
Temporally  Developing  Mixing  Layer,”  Manuscript  submitted  for 
Presentation  at  the  Twenty-First  International  Symposium  on 
Combustion,  Munich,  West  Germany,  August  1986  (Appendix  A). 
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2.  Givi,  P. ,  and  Jou,  W.-H.,  "Mixing  and  Chemical  Reactions  in  a 
Spatially  Developing  Mixing  Layer, "  Manuscript  accepted  for 
presentation  at  the  Spring  Technical  Meeting  of  the  Combustion 
Institute,  Central  States  Section,  Cleveland,  Ohio,  May  5-6,  1986 
(Appendix  B). 

INTERACTIONS 

In  addition  to  the  above-mentioned  written  manuscripts,  parts  of  our 
efforts  were  presented  (or  are  scheduled  to  be  presented)  at  the  following 
seminars  for  which  no  written  manuscripts  were  required: 

1.  Givi,  P. ,  "Direct  Numerical  Simulation  of  an  Unpremixed  Jet  Flame," 
AFOSR/ONR  Contractors  Meeting  on  Turbulent  Combustion,  California 
Institute  of  Technology,  Pasadena,  July  23-25,  1985. 

2.  Jou,  W.-H.,  "Numerical  Simulation  of  Chemically  Reacting  Flows," 
Invited  Presentation,  JANAF  Workshop  on  "Coherent  Structures  in 
Turbulent  Flows",  October  1985. 

3.  Givi,  P. ,  "Direct  Numerical  Simulations  of  a  Reacting  Mixing  Layer: 
Bringing  Fluid  Mechanics  Back  to  Turbulence  Research,"  Invited 
Seminar,  Spring  Seminar  Series,  Department  of  Mechanical  and 
Aerospace  Engineering,  State  University  of  New  York,  Buffalo, 

March  20,  1986. 
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P.  GIVI,  W.-H.  JOU  and  R.  W.  METCALFE 

Abstract 

Nonequilibrium  effects  leading  to  the  local  quenching  of  a  diffusion  flame 
have  been  investigated  by  examining  the  evolution  of  large-scale  structures  in 
a  two-dimensional  temporally  developing  mixing  layer.  Pseudospectral  calcula¬ 
tions  of  a  temperature-dependent,  nonpremixed,  reacting  shear  layer  indicate 
that  the  primary  important  parameter  to  be  considered  for  the  flame  extinction 
is  the  local  instantaneous  scalar  dissipation  rate,  conditioned  at  the  scalar 
stoichiometric  value  (Xs*).  At  locations  where  this  value  is  increased 
beyond  a  critical  value  (X^),  the  local  temperature  decreases  and  the 
instantaneous  reaction  rate  drops  to  zero.  This  is  consistent  with  the  results 
of  the  perturbation  methods  employing  large  activation  energy  asymptotics  for 
the  study  of  flame  extinction  in  nonpremixed  flames. 


1.  Introduction 

A  diffusion  flame  is  characterized  by  a  chemical  reaction  time  that  is 
usually  much  smaller  than  a  characteristic  diffusion  time.  The  chemical 
reactions  occur  in  a  narrow  zone  between  the  fuel  and  the  oxidizer,  where  the 
concentrations  of  both  reactants  are  very  small.  The  rate  at  which  the 
reactants  flow  into  the  reaction  zone,  and  therefore  the  characteristic  time 
of  the  flame,  is  dependent  on  the  hydrodynamics  of  the  particular  flow.  As 
the  characteristic  time  of  the  chemical  reactions  decreases,  this  reaction 
region  becomes  infinitesimally  thin.  In  this  limit  the  chemical  reaction  zone 
approaches  a  "flame  sheet"  (local  chemical  equilibrium*),  where  the  concen¬ 
trations  of  both  reactants  are  very  low  and  the  rate  of  combustion  is  governed 
by  the  rate  at  which  fuel  and  oxidizer  flow  into  the  reaction  zone.  The  flame 
sheet  assumption  is  justified  by  the  very  fast  chemical  reaction  rate  of  the 
diffusion  flame.  This  assumption  significantly  reduces  the  complexity  of  the 
problem  since  it  eliminates  the  analysis  associated  with  the  chemical  kinetics. 
For  many  flows  whose  characteristic  time  scale  of  chemical  reaction  is  much 


smaller  than  the  hydrodynamic  (convective-diffusive)  time  scale,  the  assumption 

of  local  chemical  equilibrium  adequately  predicts  the  location  and  the  shape 
2 

of  the  flame.  One  important  feature  of  the  calculations  based  on  the  fast 

2 

chemistry  model  is  the  introduction  of  a  passive  Shvab-Zeldovich  conserved 
scalar  variable  (Z),  which  is  independent  of  the  chemical  kinetics.  From  this 
quantity,  the  evolution  of  the  concentration  fields  of  both  the  reactants  and 
the  products  can  be  computed. 

In  turbulent  flows,  however,  the  local  characteristic  flow  time  scales 
vary  considerably.  As  a  result,  many  important  and  interesting  problems  that 

cannot  be  analyzed  by  local  chemical  equilibrium  assumptions  are  introduced. 

.  3 

Experimental  studies  of  Tsuji  show  that  as  the  local  characteristic 

diffusion  time  becomes  shorter  and  approaches  the  order  of  magnitude  of  the 

chemical  time  scale,  the  details  of  the  chemical  reactions  cannot  be  neglected. 

If  the  flow  of  reactants  into  the  reaction  zone  increases  further,  causing  the 

diffusion  time  scale  to  be  reduced  more,  the  chemical  reaction  will  not  be 

able  to  keep  pace  with  the  further  supply  of  reactants.  The  reaction  rate 

4 

will  be  reduced,  and  local  quenching  occurs.  As  shown  by  Peters,  further 
reduction  of  the  diffusion  time  scale  leads  to  lift-off  and  the  blow-off  of 
the  entire  flame.  As  noted  by  both  Tsuji  and  Peters,  the  consideration  of 
flame  extinction  cannot  be  explained  by  the  flame  sheet  model,  which  assumes 
an  infinitely  fast  chemical  reaction  rate.  Therefore,  in  order  to  address  the 
quenching  phenomenon,  the  structure  of  a  finite  reaction  rate  zone  must  be 
studied . 

There  are  a  number  of  possible  ways  to  account  for  the  influence  of  finite 
rate  chemistry  on  a  diffusion  flame.  For  laminar  flows,  a  full  set  of 
conservation  equations  with  one-  or  two-step  kinetics  has  been  solved  by 
numerical  methods.^  With  simplified  kinetics,  perturbation  methods  based  on 


large  Damkohler  number  asymptotics^*^  or  on  large  activation  energy 
8  9 

asymptotics  *  have  proven  useful.  Damkohler  number  asymptotics  give  some 
estimate  of  the  broadening  of  the  reaction  zone  in  nonequilibrium  situations. 
The  analysis  of  extinction  and  ignition  cannot  be  addressed  by  this  method, 
however,  but  can  be  analyzed  by  activation  energy  asymptotics.  (For  a 
detailed  summary  and  review  of  the  different  perturbation  techniques,  see 
Ref.  2.) 

9 

Linan  has  employed  a  method  of  matched  asymptotic  expansions  in  the 

limit  of  large  activation  energy  in  an  attempt  to  describe  the  interaction 

between  the  hydrodynamics  and  chemical  reactions  in  the  reaction  zone  of  a 

counter flowing  laminar  diffusion  flame.  It  has  been  shown  that  activation 

energy  asymptotics  are  very  useful  in  predicting  flame  ignition  and  extinction 

characteristics  in  such  flows.  This  technique  was  extended  to  turbulent  flows 
4  10 

by  Peter s  ’  by  considering  turbulent  diffusion  flames  as  an  ensemble  of 

laminar  diffusion  flamelets.11  By  introducing  a  local  coordinate  system 

that  moves  with  the  stoichiometric  flame  sheet,  and  employing  a  perturbation 

4 

of  the  large  activation  energy  asymptote,  Peters  was  able  to  recognize  that 
the  primary  "nonequilibrium"  parameter  for  the  analysis  of  the  flame 
extinction  is  the  dissipation  rate  of  the  scalar  quantity  evaluated  at 
stoichiometric  conditions  (Xgt).  This  quantity  is  viewed  as  the  inverse  of 
the  diffusion  time  scale.  If  this  parameter  is  increased  beyond  a  critical 
limit  (X  ),  the  heat  conducted  to  both  sides  of  the  diffusion  flamelets 

q 

cannot  be  balanced  by  the  heat  production  due  to  chemical  reaction.  At  the 

critical  value  of  the  dissipation,  the  finite  rate  kinetics  balance  the 

3  12 

diffusion.  Some  numerical  calculations,  performed  by  Liew  et  al.  and 

.  .  13 

compared  with  experimental  data,  show  also  that  as  the  maximum  value  of  the 


dissipation  (X  )  increases,  the  value  of  the  maximum  temperature  decreases, 
max 

the  reaction  eventually  ceases,  and  the  flame  is  locally  quenched. 

These  results  indicate  that  the  dissipation  rate  of  the  conserved  scalar 
is  a  useful  characteristic  to  study  in  the  analysis  of  nonequilibrium  effects 
leading  to  flame  extinction.  In  turbulent  flows,  however,  the  quantity  X 

is  a  strongly  fluctuating  quantity  that  has  not  yet  been  numerically  computed. 

14 

Instead,  statistical  approaches  have  been  chosen  by  Peters  and  Williams 
and  Janicka  and  Peters^  employing  different  turbulence  closures  in  pre¬ 
dicting  the  lift-off  height  of  a  round  diffusion  jet  flame.  Results  obtained 
using  these  methods  were  compared  with  experimental  data^  for  both  methane 
and  natural  gas  flames  and  show  some  correct  order-of-magnitude  predictions. 
These  results  are  encouraging;  however,  with  the  availablity  of  larger 
computers,  it  is  now  possible  to  employ  a  more  accurate  treatment  of  the  flame 
extinction  and  local  flame  quenching  that  occur  in  nonpremixed  flames.  It 
would  be  very  useful  to  correlate  the  flame  extinction  with  the  nonequilibrium 
parameter  Xgt.  Also,  it  would  be  very  interesting  to  look  at  the  structure 
of  the  diffusion  flames  at  the  point  of  extinction.  These  are  the  objectives 
of  this  paper. 

The  present  paper  employs  a  direct  numerical  simulation  approach  to 
investigate  the  problem  of  local  flame  extinction  in  a  time-dependent,  two- 
dimensional  mixing  layer.  In  this  flow,  the  governing  equations  are  solved  by 
an  accurate  numerical  method  without  a  closure  model.  The  time-dependent  flow 
field  can  be  analyzed  statistically  to  understand  the  underlying  physics  much 
as  an  experimentalist  does  with  laboratory  data.  The  direct  numerical 
simulation  technique  has  recently  been  successfully  applied  to  chemically 
reacting  flows.  Riley  et  al.^  considered  the  three-dimensional  temporally 
growing  mixing  layer  under  the  simplest  possible  assumption  of  a  constant 


rate,  non-heat-releasing  chemical  reaction.  The  main  contribution  of  the  work 

is  the  understanding  of  the  effects  of  three-dimensional  mixing  and  the  diffu- 

18 

sion  of  the  species  on  the  chemical  reactions.  McMurtry  et  al.  considered 
the  effects  of  the  chemical  heat  release  and  the  resulting  density  variation 
on  the  fluid  motion  for  a  two-dimensional  mixing  layer.  The  fluid  dynamics 
and  the  chemical  reaction  are  truly  coupled  in  this  work,  and  the  interplay 
between  the  two  are  discussed.  However,  the  assumption  of  a  constant  chemical 
reaction  rate  is  still  employed.  In  the  present  work,  we  intend  to  understand 
the  flame  extinction  problem  through  a  two-dimensional  simulation  of  a  mixing 
layer.  In  particular,  the  role  of  large-scale  features  of  the  turbulent  flow 
in  flame  extinction  is  studied.  Due  to  limitations  of  numerical  accuracy, 
only  moderate  Reynolds  and  Damkohler  numbers  are  computed.  For  flame  extinc¬ 
tion  in  a  three-dimensional  turbulent  flow,  three-dimensional  simulations  will 
be  performed  in  the  future. 

In  the  following  section,  the  geometrical  configuration  of  the  problem  is 
given,  together  with  the  formulation  of  the  problem.  In  Section  3,  some 
sample  results  are  presented.  Finally,  in  Section  4,  conclusions  are 
presented . 

2.  Problem  Formulation 

For  simplicity  of  numerical  calculations,  we  consider  a  two-dimensional 
shear  layer  as  shown  in  Fig.  1.  The  reactant  flows  to  the  right  on  the 
upper  stream,  and  the  other  reactant,  CA,  flows  to  the  left  on  the  lower 
stream.  The  flow  is  considered  periodic  in  the  horizontal  direction  Cx). 
Although  the  splitter  plate  flow  evolves  spatially  downstream  and  the  numerical 
simulations  evolve  temporally ,  important  similarities  in  the  dynamics  of  these 
two  flows  make  it  useful  to  study  accurate  numerical  simulations  of  the 


temporally  growing  mixing  layers.  According  to  the  model,  a  flow  quantity 
averaged  in  the  x-direction  can  be  related  to  the  time  average  of  the  6ame 
quantity  at  a  fixed  location  in  a  splitter  plate  configuration.  These  average 
quantities  are  dependent  on  the  transverse  coordinate  (y)  and  the  time  (t). 
Again,  the  inverse  Galilean  transformation  relates  the  time  t  to  the  stream- 
wise  location  in  a  splitter  configuration.  This  allows  us  to  use  accurate 
pseudospectral  numerical  methods;  these  methods  are  discussed  in  Ref.  17  and 
will  not  be  repeated  here. 

Initially,  the  two  chemical  species  A  and  B  are  not  premixed.  The 
chemical  reaction  between  the  two  species  is  assumed  to  be  single-step  and 
irreversible  and  to  obey  the  temperature-dependent  Arrhenius  law.  It  is 
assumed  that  the  heat  release  rate  is  low  so  that  the  effects  of  the  chemical 
heat  release  on  the  flow  field  can  be  neglected.  This,  together  with  the 
further  assumptions  of  a  low  Mach  number  flow,  result  in  a  constant-density 
flow  formulation. 

The  velocity  field  (U)  satisfies  the  incompressible  Navier-Stokes 
equations  normalized  by  the  density,  which  is  taken  to  be  unity: 

7*U  =  0  (1) 

3U  , 

l(u,  v)  »  +  u-yu  -  vru  =  -  vp  (2) 

The  species  concentration  fields  are  governed  by 

L(Ca,  D)  =  L(Cb,  D)  =  -  W  (3) 

In  this  equation,  the  stoichiometric  coefficient  is  set  equal  to  unity,  and 
the  reaction  rate  follows  the  Arrhenius  law: 

W  =  Af  e‘E/RT  CACB  (4) 


.v.v.v 


•  ^  .  •7'-' .  -  .  j  i.  -Vv.1. 


»  «  *  -  . 


V 
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where  A^  is  the  preexponential  frequency  factor,  E  is  the  activation  energy,  R 
is  the  universal  gas  constant,  and  T  is  the  temperature,  which  is  governed  by 


its  own  transport  equation. 


C  L(T,  K)  =  QW 
v 


To  identify  the  nondimen sional  groupings,  the  variables  are  normalized  by  a 
velocity  scale  (mean  velocity  difference  across  the  layer  &J),  a  temperature 
scale  (free-stream  T^) ,  concentration  scales  (free-stream  CA<X)  ®nd  CB<B) ,  and  a 
length  scale  (X  *  2.249  C,  where  O  is  the  distance  from  the  plane  of  symmetry 
to  the  transverse  plane,  where  the  mean  velocity  rises  to  half  of  its  free- 
stream  value).  The  value  of  2.249  was  chosen  so  that  the  size  of  the 
computational  domain  would  be  equal  to  the  wavelength  of  the  most  unstable 
mode  and  its  normalized  value  would  be  an  integer  multiple  of  2lt. 
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Assuming  equal  free-stream  concentrations,  i.e.,  CAoo  *  CBeo,  the  nondimen- 
sionalized  equations  become 


V«U*  *  0 


*,  *  „  *  _  *  i  * 

L  (U  ,  Re)  *  +  U  *7U  -  jj^-  7  U  = 

*  *  *  *  * 

L  (CA,  Re  Sc)  -  L  (C„,  Re  Sc)  -  -  W 
A  D 

*.  ★  .  * 

L  (Cp,  Re  Sc)  -  W 

L* ( T* ,  Re  Pr)  -  Ce  W* 
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where 
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Ze 
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-Ze/T  *  * 
Da  e  '  CAC„ 
A  B 


AfC» 

AU 

r 


Damkohler  number 


E 

RT 

o 

&JX 

V 

_2_ 

C  T 


Zeldovich  number 

Reynolds  number 
Heat  release  parameter 

Schmidt  number 
Prandtl  number 


Normalization  of  the  equations  yields  the  important  nondimens ional  groupings 
that  appear  in  the  formulation:  the  Reynolds  number,  the  Damkohler  number, 
the  Zeldovich  number,  the  heat  release  parameter,  the  Schmidt  number,  and  the 
Prandtl  number.  The  values  of  Sc  and  Pr  can  be  set  equal  to  unity,  since  this 
is  approximately  correct  for  gaseous  diffusion  flames.  This  results  in  a 
Lewis  number  of  unity.  The  value  of  other  nondimensional  parameters  is 
limited  by  the  available  resolution  of  the  numerics  employed  in  the 
computations . 


Shvab-Zeldovich  Variable 

For  the  two-feed  diffusion  flame  considered  here,  it  is  possible  to 

consider  only  three  scalar  quantities  rather  than  solving  for  all  four  scalar 

variables.  The  Shvab-Zeldovich  variable  (Z)  and  the  product  concentrations 

19 

(G),  following  Givi  et  al.,  are  defined  as 

*  * 

2  ■  CA  *  CP 
* 


G«Cp 


-  >N> 


•  ■  \.V\‘A-yyv-.v.’’s 


The  transport  equations  governing  Z,  G,  and  T  follow: 


*/  » 

L  (Z)  =  0 

★ 

L*(G>  -  Da  e”Ze'T  (1  -  Z  -  G)  (Z  -  G)  (10 

* 

L*(T*)  *  Ce  Da  e~Ze/T  (1  -  Z  -  G)  (Z  -  G) 

One  could  further  relate  the  normalized  temperature  and  product  concentrations 

to  the  normalized  Shvab-Zeldovich  variable,  eliminating  the  need  for  the  solu- 

* 

tion  of  the  equation  governing  G,  if  the  initial  conditions  for  T  and  G  could 
be  related.*  For  the  moderate  Damkohler  numbers  employed  in  these  calculations 
however,  the  reaction  zone  at  the  interface  of  the  two  reactants  was  numeri¬ 
cally  ignited  by  increasing  the  temperature  around  the  reaction  surface, 
resulting  in  different  initial  normalized  profiles  for  the  product  and  tempera¬ 
ture.  Therefore,  for  the  treatment  of  the  nonpremixed  system  we  plan  to 

investigate,  the  full  solution  of  hydrodynamic  variables  (U  ,  V  ,  P  ),  scalar 

★ 

variables  (Z,  G)  and  temperature  (T  )  is  required. 


Initial  and  Boundary  Conditions 

The  initial  and  boundary  conditions  for  the  velocity  field  are  given 
20 

elsewhere  and  are  only  summarized  here.  In  terms  of  the  stream  function, 
the  initial  condition  is  given  by 


* 


♦ 

mean 


+  i|> 


sh 


(11 


where  is  the  stream  function  associated  with  a  mean  hyperbolic 

tangent  velocity  profile,  is  the  stream  function  for  the  most  unstable 

mode  of  this  mean  velocity  profile,  and  is  the  stream  function  for  the 

first  subharmonic  of  the  most  unstable  mode.  The  properties  of  these  modes 

21 

have  been  evaluated  from  linear  stability  theory.  The  fundamental  mode  in 


this  mixing  layer  produces  a  single  vortex  rollup 
added  in,  a  second  rollup,  or  pairing,  can  occur. 


3.  Presentation  of  Results 


■  v.v.v.v  ,-v1-  r. 


A  pseudospectral  numerical  code  developed  by  McMurtry  et  al .  was  modi¬ 
fied  for  the  calculations  of  the  incompressible  reacting  flow  considered  here. 
Computations  were  performed  in  a  square  domain  with  the  size  0  <  x  <  2ltX, 

-irX  <  y  <  hX  for  the  single  vortex  rollup  and  the  domain  0  <  x  <  4ttX, 

-2irX  <  y  <  2irX  for  the  double  rollup  computations.  The  spatial  resolution  is 
64x64  Fourier  modes.  The  values  of  the  Reynolds  and  Damkohler  numbers  were  set 
equal  to  200  and  10,  respectively,  so  that  the  simulation  would  be  accurately 
resolved  on  the  64x64-point  grid  employed  here.  The  value  of  the  Zeldovich 
number  based  on  the  free-stream  temperature  was  set  equal  to  20.  The  value  of 
the  heat  release  parameter  was  selected  to  be  small  enough  so  that,  when  multi¬ 
plied  by  the  Damkohler  number,  it  would  be  resolvable  by  the  numerics  and  so 
that  it  would  also  be  reasonable  to  neglect  the  effects  of  heat  release  on  the 
flow  field.  However,  this  value  should  be  large  enough  to  allow  the  effects 

of  temperature  dependence  on  the  chemical  reaction  term.  A  value  of  Ce  *  8  is 

★ 

satisfactory,  giving  a  maximum  (adiabatic)  temperature  of  Tacjiabatic  = 

Because  of  diffusion,  the  actual  maximum  temperature  at  the  reaction  surface 
will  be  less  than  the  above  flame  temperature  obtained  under  adiabatic 
conditions . 

The  velocity  profile  for  incompressible  periodic  mixing  layers  has  already 
been  documented*^  and  will  not  be  discussed  here.  Instead,  the  chemical 
quantifier  that  are  of  interest  in  flame  extinction  problems  will  be  presented. 
In  this  section,  the  instantaneous  profiles  of  the  concentration  fields,  the 
temperature  field,  and  the  local  reaction  rates  are  presented  and  discussed  in 
detail. 

The  profiles  of  the  conserved  Shvab-Zeldovich  scalar  variable  are  chosen 
for  the  purpose  of  flow  visualization.  In  Figs.  2  and  3,  we  present  the  time 


11 


sequence  development  of  this  profile  for  Cases  1  and  II,  respectively. 

Initially,  the  perturbation  associated  with  the  fundamental  mode  grows  until  a 
* 

time  of  t  =9,  where  the  first  vortex  rollup  occurs.  This  results  in  a 
very  steep  gradient  of  the  Shvab-Zeldovich  variable  near  the  braids  of  the 
vortex  rollup  due  to  the  high  strain  field  there.  Proceeding  further  in  time 
results  in  diffusion  of  the  core  of  the  vortex  with  no  additional  rollup. 

Adding  the  subharmonic  associated  with  this  unstable  mode  results  in  a  second 

.  ...  .  * 
rollup  (Fig.  3)  that  initiates  at  a  time  of  about  t  =12  and  is  completed 

* 

by  t  =  24.  As  shown  in  these  figures,  the  effect  of  the  vortex  dynamics  is 

to  increase  the  strain  rate  at  the  braids  of  the  rollup. 

If  a  fast  chemistry  model  was  assumed  to  describe  the  chemical  reactions, 

the  surface  of  Z  =  Z  „  =  0.5  would  correspond  to  the  flame  location.  The 

st  r 

vortex  dynamics  plays  an  important  role  in  increasing  the  area  of  the  flame 
surface  (Zgt  surface)  and  in  increasing  the  total  amount  of  reaction  product 
compared  to  the  case  where  there  is  no  rollup  and  the  chemical  reaction  is 
only  limited  to  the  interdiffusion  of  the  scalars  at  the  reaction  surface.  In 
the  present  calculations,  the  increase  of  the  equilibrium  flame  sheet  surface 
due  to  rollup  (Case  1,  Fig.  2c)  is  153  percent  in  comparison  with  the  unforced 
case.  Merging  of  the  vortices  (Case  II,  Fig.  3c)  increases  this  surface  by 
300  percent  when  compared  with  the  unforced  case. 

The  time  sequence  of  the  product  concentration  is  also  shown  on  Figs.  4 
and  5  for  Cases  1  and  II,  respectively.  The  mechanism  of  vortex  rollup  is  to 
draw  the  reactants  from  the  two  streams  to  the  reaction  surface.  The  amount 
of  products,  along  the  Z  =  Zgt  surface,  reaches  a  maximum  at  the  core  of  the 
vortex  and  a  minimum  at  the  braids.  The  same  trends  are  shown  by  the  tempera¬ 
ture  profiles  (Figs.  6  and  7),  where  the  maximum  temperature  occurs  at  the 
core  (the  hottest  location),  where  vorticity  is  also  highest,  and  decreases  in 


the  braids,  where  the  gradients  of  the  scalar  are  at  their  higher  values.  At 

the  time  where  the  computation  is  stopped,  the  maximum  product  concentration 

has  not  yet  reached  its  equilibrium  value.  There  still  remain  reactants,  even 

at  the  core  of  the  mixing  layer,  and  the  flame  shortening  phenomenon,  which  is 

usually  caused  by  the  depletion  of  the  reactants,  has  not  been  yet  observed. 

The  reason  for  this  behavior  of  the  chemical  product  and  temperature  can  be 

explained  by  the  contour  plots  of  the  instantaneous  reaction  rate  W,  presented 

in  Fig.  8  for  Case  I.  As  shown  on  this  figure,  the  reaction  rate  initially  is 

very  much  uniform  along  the  species  interdiffusion  zone  and  is  maximized  at 

the  reaction  surface,  where  the  reactants  are  in  contact.  At  later  times,  the 

reaction  rate  decreases,  and  at  points  where  the  strain  rate  is  sufficiently 

large,  the  reaction  rate  goes  to  zero  and  the  flame  locally  quenches.  This 

mechanism  of  extinction  of  the  diffusion  flame  is  consistent  with  the  experi- 

3 

mental  observations  of  Tsuji.  In  the  braids,  as  the  inverse  of  the 

diffusion  time  scale  (1/X  )  decreases  (as  the  result  of  the  vortex  rollup), 

s  t 

the  supply  of  the  reactants  is  greater  than  the  chemical  reaction  can  keep 
pace  with.  Notice  from  Figs.  2  through  5  that  the  reactants  are  well  mixed  on 
the  braids,  but  the  reaction  rates  at  these  places  are  zero.  This  nonequili¬ 
brium  phenomenon  can  be  explained  by  examining  the  reaction  rate  equation 
[Eq.  (4)]. 

Unlike  the  temperature-independent  reacting  mixing  layer  calculation 
.  22 

reported  by  Givi  et  al.,  the  maximum  value  of  the  reaction  rate  does  not 
necessarily  correspond  to  the  maximum  value  of  the  product  of  the  reactants 
concentration,  since  the  effects  of  the  temperature  variations  influence  the 
local  conversion  rate  of  the  chemistry.  If  the  flame  temperature  goes  below  a 
critical  characteristic  temperature,  the  flame  becomes  very  rich  with  both 
reactants.  However,  since  the  value  of  the  dissipation  rate  is  greater  than 


the  critical  value,  this  premixed  region  of  the  reactants  cannot  be  reignited 
from  the  high  temperature  zone  of  the  -core.  A  higher  dissipation  rate  results 
in  a  further  decrease  in  temperature  until  it  becomes  equal  to  the  background 
temperature.  At  this  point,  the  vortex  has  reached  the  boundary  of  the 

computational  domain  and,  therefore,  further  computation  is  not  realistic. 

•  •  A 

As  mentioned  earlier,  Peters  developed  a  theory  based  on  the  methodology 
9 

of  Lxnan  to  obtain  a  criterion  for  the  critical  value  of  the  dissipation 
rate  for  an  irreversible  one-step  chemical  reaction.  Employing  large  activa¬ 
tion  asymptotics,  he  derived  an  equation  for  X^  that  depends  on  the  frequency 
factor  of  the  kinetics  term,  the  activation  energy,  Zgt,  and  the  temperature 
at  the  stoichiometric  value.  A  direct  quantitative  comparison  of  X^,  obtained 
in  our  numerical  simulations  (X^  is  approximated  here  as  the  value  of  the 

instantaneous  dissipation  rate  evaluated  at  the  stoichiometric  surface  Z  ) 

st 

with  that  obtained  under  the  large  activation  energy  asymptotics  is  not 
expected  to  agree  exactly  for  several  reasons.  First,  in  asymptotic  analysis, 
the  flame  thickness  is  very  small  and  is  only  broadened  near  the  reaction  zone 
by  a  small  parameter.  In  the  direct  numerical  simulations  reported  here,  the 
reactants  have  a  fairly  "thick"  overlap  region  in  the  reaction  zone  due  to  the 
numerical  accuracy  limitations.  Furthermore,  the  initial  conditions  for  the 
temperature  profile  employed  here  near  the  reaction  zone  are  not  the  same  as 
the  uniform  initial  temperature  distributions  of  Peters.  The  value  of  the 
initial  temperature  is  very  important  in  the  calculation  of  the  adiabatic  flame 


temperature  employed  in  the  asymptotic  analysis.  This  was  also  discussed  by 
12 

Liew  et  al.  in  their  numerical  prediction  of  the  flame  extinction.  However, 


making  the  assumption  that  the  initial  concentration  of  the  reactant  and  the 
temperature  of  the  reactants  are  at  some  average  values  in  the  range  between 
the  free  stream  and  the  reaction  zone,  correct  order-of-magnitude  asymptotic 


results  are  verified  by  the  results  of  the  simulations.  The  temperature  pro- 

* 

files  presented  in  Fig.  6  show  that  by  a  time  of  t  =  15,  the  value  of  the 
temperature  at  the  braids  falls  to  about  one-fifth  the  flame  temperature 
(approximately  equal  to  the  free-stream  temperature).  This  corresponds  to 

* 

local  quenching.  The  corresponding  computed  normalized  dissipation  rate  X 
“k 

(Xgt  =  XgtA/AU)  at  the  point  of  quenching,  for  the  given  kinetics  data,  is 
about  7  for  both  Cases  I  and  II.  The  estimated  corresponding  value  for  the 
dissipation  rate,  computed  with  the  asymptotic  methods  described  above,  is 
about  8,  which  is  in  good  agreement  with  the  results  of  the  numerical 
simulation. 


4.  Conclusions 

A  pseudospectral  algorithm  has  been  used  for  the  numerical  calculations  of 
a  chemically  reacting,  temperature-dependent  mixing  layer.  The  nonequilibrium 
effects  leading  to  the  local  flame  quenching  have  been  simulated  in  a  case 
with  fairly  large  Zeldovich  number  and  moderate  Reynolds  and  Damkohler  numbers. 
It  has  been  found  that  the  rollup  of  an  unstable  shear  layer  creates  regions 

with  high  dissipation  rates  at  the  braids,  where  local  flame  extinction  occurs 

4  14 

according  to  the  theory  of  Peters  ’  .  The  temperature  contour  shows  that 

the  temperature  drops  to  a  value  close  to  ambient  at  the  braids  and  the  reac¬ 
tion  rate  reduces  to  zero,  even  though  the  reactants  are  well  mixed  there. 
Within  the  limitations  of  numerical  resolution,  comparison  of  the  critical 
scalar  dissipation  rate  with  that  obtained  by  large  activation  energy 
asymptotics  shows  reasonable  agreement.  We  are  presently  expanding  this  work 
to  investigate  a  spatially  developing  mixing  layer  and  three-dimensional 
turbulent  flows.  This  should  further  illuminate  the  basic  mechanisms  of  such 


phenomena  as  the  diffusion  jet  flame  lift-off 


Preexponential  factor 
Concentration 


Nomenclature 


Specific  heat 

Heat  release  parameter 

Molecular  diffusivity 

Damkohler  number 

Activation  energy 

Normalized  product  concentration 

Thermal  conductivity 

Convective-diffusive  operator 

Pressure 

Prandtl  number 

Heat  of  reaction 

Universal  gas  constant 

Reynolds  number 

Schmidt  number 

Temperature 

Time 

Streamwise  velocity  component 
Cross-stream  velocity  component 
Reaction  rate 
Dissipation 
Streamwise  coordinate 
Cross-stream  coordinate 
Shvab-Zeldovich  variable 


Zeldovich  number 
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Figure  Captions 


FIG  1.  Problem  geometry. 

FIG  2.  Plots  of  Shvab-Zeldovich  variable  contours  for  a  sequence  of 
times  (Case  1). 

FIG  3.  Plots  of  Shvab-Zeldovich  variable  contours  for  a  sequence  of  times 
(Case  11). 

FIG  4.  Plots  of  normalized  product  concentration  contours  for  a 
sequence  of. times  (Case  1). 

FIG  5.  Plots  of  normalized  product  concentration  contours  for  a 
sequence  of  times  (Case  II). 

FIG  6.  Plots  of  normalized  temperature  contours  for  a  sequence 
of  times  (Case  I). 

FIG  7.  Plots  of  normalized  temperature  contours  for  a  sequence 
of  times  (Case  II). 

FIG  8.  Plots  of  normalized  reaction  rate  contours  for  a  sequence 
of  times  (Case  1). 
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PIG  7.  Plots  of  normalised  temperature  contours  for  a  sequence 
of  times  (Case  II). 
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FIG  8.  Plots  of  normalized  reaction  rate  contours  for  a  sequence 
of  times  (Case  I). 
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ABSTRACT 

Direct  numerical  simulations  have  been  employed  to  study  the  effects  of 
large  coherent  structures  in  a  perturbed,  time-dependent,  spatially  developing, 
reacting  mixing  layer.  It  is  shown  that  the  vortex  dynamics  has  a  significant 
influence  on  the  enhancement  of  the  chemical  reactions,  convolution  of  the 
reaction  surface,  and  increase  of  product  formation  in  the  reaction  zone  of  the 
layer.  Furthermore,  examination  of  some  of  the  statistical  quantities  obtained 
from  the  results  of  direct  numerical  simulations  indicates  the  asymmetric 
nature  of  the  mechanism  of  mixing  processes  in  the  mixing  zone  of  the  layer. 

1.  INTRODUCTION 

The  advancement  of  supercomputer  technology  in  recent  years  has  had  a  major 
influence  on  Increasing  our  understanding  of  the  mechanisms  of  many  complex 
physical  phenomena  such  as  turbulence.  The  availability  of  high-speed,  large- 
memory  computers  has  made  it  possible,  in  many  cases,  to  simulate  turbulent 
flow  problems  directly  by  solving  the  appropriate  basic  governing  transport 
equations  without  any  need  for  additional  "turbulence  modeling. "  Such  modeling 
is  required  when  these  equations  are  statistically  averaged.  In  a  complex 
system  such  as  a  chemically  reacting  flow,  modeling  is  extremely  difficult 
because  of  our  lack  of  knowledge  on  the  detailed  dynamics  of  the  flow.  Without 
using  turbulence  modeling,  direct  numerical  simulations  can  provide  the  de¬ 
tailed  Information  needed  to  Increase  our  understanding  of  the  physical  pro¬ 
cesses  Involved  in  turbulent  flows.  An  increase  in  understanding  various 
aspects  of  the  physics  of  turbulent  reacting  flows  is  evident  in  recent  works 
(e.g.,  Riley  et  al.,  1986). 


Direct  numerical  simulations  o£  turbulent  reacting  flows  are  presently 
applied  to  flows  with  limited  dynamical  and  chemical  parameter  ranges. 
Continuing  efforts  in  evaluating  the  advantages  and  also  the  limitations  of 
this  technique  are  required  before  the  knowledge  acquired  can  be  applied  to 
flows  with  practical  parameter  ranges.  The  application  of  direct  numerical 
simulations  to  a  flow  with  simplified  geometrical  configurations  and  well** 
controlled  conditions  is  certainly  the  first  step  in  understanding  the  physics 
and  in  evaluating  the  capabilities  of  the  method. 

A  mixing  layer,  in  which  two  parallel  streams  with  different  velocities 
begin  to  mix  and  react  downstream  of  the  trailing  edge  of  a  splitter  plate 
partition,  presents  an  ideal  configuration  for  studying  the  physical  processes 
in  the  mixing  and  reaction  zones  that  exist  in  real  diffusion-controlled  com¬ 
bustors.  Direct  numerical  simulations  have  been  applied  in  studying  a  model 
problem  of  such  flows  with  encouraging  success  (McMurtry  et  al.,  1985;  Riley 
et  al.,  1986).  However,  most  of  the  previous  work  has  employed  a  temporally 
developing  mixing  layer  as  a  model.  Temporally  growing  mixing  layers  refer  to 
flows  that,  in  a  Galilean-transformed  frame  of  reference,  are  assumed  to  evolve 
in  time  rather  than  in  space.  This  is  a  good  approximation  for  flows  in  which 
the  difference  between  the  freestream  velocities  of  two  streams  is  much  less 
than  the  averaged  velocity.  This  approximation  allows  simplifications  in  the 
formulation  of  periodic  boundary  conditions.  A  pseudospectral  numerical  method 
using  Fourier  series  can  then  be  employed  to  solve  the  appropriate  transport 
equations. 

Comparison  of  the  simulation  results  obtained  for  a  temporally  developing 
mixing  layer  with  experimental  results  is  encouraging.  The  comparison  shows 
that  there  is  qualitative  agreement  between  numerical  simulations  and  labor¬ 
atory  data  for  the  average  reactant  concentrations  and  the  concentration 
correlations.  However,  there  are  some  basic  phenomena  that  cannot  be  explained 
by  temporal  simulations,  such  as  (1)  asymmetric  mixing  observed  in  mixing  layer 
experiments  (Koochesfahani,  1963),  and  (2)  stablization,  quenching,  lift-off, 
and  blowout  of  turbulent  jet  flames  (Peters  and  Williams,  1983).  A  fundamental 
understanding  of  these  phenomena  requires  detailed  calculations  of  spatially 
growing  flows.  Here,  we  have  developed  a  numerical  technique  that  is  capable 
of  simulating  a  spatially  developing  mixing  layer  as  a  first  step  to  further 
develop  the  technique  of  direct  numerical  simulations. 
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This  paper  presents  some  preliminary  results  of  our  two-dimensional  flow 
simulations  using  the  resulting  code  in  an  attempt  to  address  the  first  of  the 
above-mentioned  phenomena  associated  with  spatially  growing  flows.  The  flow 
fields  under  consideration  are  dominated  by  large-scale  fluctuations  associated 
with  coherent  structures.  The  chemical  reaction  that  occurs  between  the  reac¬ 
tants  on  two  sides  of  the  layer  is  assumed  to  be  infinitely  fast,  so  that 
complications  due  to  complex  kinetics  need  not  be  considered.  Therefore,  the 
transport  of  a  conserved  Shvab-Zeldovich  scalar  variable  is  computed.  There 
is  also  a  nonparallel  growth  mechanism  of  the  layer  mainly  due  to  diffusion, 
vortex  rollup,  and  subsequent  pairing  of  the  neighboring  vortices.  The  effects 
of  this  nonparallel  growth  as  well  as  the  asymmetric  mixing  in  the  shear  layer 
are  studied.  Analysis  of  flame  extinction  and  lift-off  phenomena  requires 
consideration  of  nonequllibrium  chemistry.  This,  as  well  as  the  effects  of  a 
three-dimensional  random  turbulence  field,  are  presently  under  study  and  will 
be  reported  in  the  future. 

In  the  next  section,  the  governing  equations,  boundary  conditions,  and 
numerical  algorithms  for  the  present  problem  are  discussed.  Some  sample 
results  are  presented  in  Section  3,  and  conclusions  drawn  based  on  these 
results  are  given  in  Section  4. 

2.  PROBLEM  FORMULATION 

The  mixing  layer  under  consideration  is  shown  in  Figure  1.  The  reactants 
are  Introduced  as  dilute  traces  in  a  carrier  passive  gas.  The  chemical  heat 
release  is  assumed  negligible.  Therefore,  the  density  of  the  mixture  can  be 
assumed  constant  and  equal  to  that  of  the  carrier  gas.  Under  these  assump¬ 
tions,  the  hydrodynamics  and  the  transport  of  chemical  species  are  decoupled. 

The  latter  is  governed  by  passive  transport  equations  with  a  given  hydrodynamic 
velocity  field. 

The  velocity  field  satisfies  the  incompressible  Navier-Stokes  equations: 

p  L(U,  v)  -  -  VP  Cl) 

V«U  -  0  (2) 

where  P  is  the  pressure,  p  is  the  density,  which  is  assumed  to  be  constant,  v 
is  the  kinematic  viscosity,  U  is  the  velocity  vector  and  L  is  the  convective 
diffusive  operator  defined  by 

/  3U\ 


Species  A  is  injected  from  the  high-speed  stream  on  the  upper  side  of  the 
plate,  and  species  B  is  Introduced  from  the  low-speed  stream.  The  chemical 
reaction  between  the  species  field  is  assumed  to  be  fast  and  irreversible, 
such  as 

A  +  B+  C  +  D  .  (4) 

The  molar  concentrations  of  the  reactants  (species  A  and  B)  and  the  product 
(say,  species  C)  satisfy  the  following  reaction-convection-diffusion  equations: 


l(ca,  D) 
L(CC,  D) 


l(cb,  D) 


-  W 


(5) 

(6) 


where  W  is  the  reaction  rate,  and  the  stoichiometric  coefficient  is  assumed  to 
be  equal  to  1.  The  molecular  diffuslvities,  D,  of  the  species  are  assumed  to 


be  equal,  i.e.,  D 


ij  ji 


D.  If  we  define  a  Shvab-Zeldovich  variable  by 


CA  +  CC 


(7) 


then  the  transport  equation  for  f  satisfies 

L(f ,  D)  -  0  .  (8) 

As  shown  by  Williams  (1985)  and  Bilger  (1980),  assuming  an  infinitely  fast 
chemical  reaction  between  species  A  and  B,  one  would  be  able  to  determine  the 
instantaneous  fields  of  all  the  species  involved. 

In  order  to  solve  the  set  of  Equations  (1)  through  (8),  boundary  and 
initial  conditions  are  needed.  For  the  cross-stream  direction,  we  assume  the 
flow  to  be  periodic  with  free-slip  boundary  conditions.  At  the  inflow,  the 
initial  conditions  for  the  Shvab-Zeldovich  variable,  f,  are  given  by  the 
following  functional  form: 


The  Initial  conditions  for  the  mean  streamwise  velocity  are  given  by  a  hyper¬ 
bolic  tangent  velocity  profile,  and  the  mean  cross-stream  velocity  component 
is  set  equal  to  zero: 

U  -  [U(Y) ,  0]  (10) 

U(Y)  -  Ux  [tanh  (2Y/A)  +  2] 


(11) 


where  Is  the  slow  stream  velocity,  and  A  la  the  vortlcity  thickness  defined 
by 

A  -  AU/(dU/dY)  .  (12 

max 

The  term  AU  is  the  difference  between  the  velocities  in  the  two  layers. 

It  Is  well-known  that  a  shear  layer  subject  to  white  noise  will  develop 
coherent  structures  with  predominantly  the  most  unstable  mode  of  its  insta¬ 
bility  waves,  and  the  subsequent  pairing  is  caused  by  the  subharmonic  distur¬ 
bance  of  the  predominant  mode.  Therefore,  in  order  to  initiate  the  vortlcity 
rollup  and  pairing,  a  small  perturbation,  in  the  form  of  the  most  unstable 
mode  and  its  subharmonics,  calculated  from  linear  stability  theory  (Michalke, 
1965;  Monkewltz  and  Huerre,  1982)  is  added  to  this  mean  profile.  In  terms  of 
stream  function,  we  have 

♦  *  +  h  *1  *  e2  *2  (i: 


where  4*  is  the  stream  function  associated  with  the  mean  flow,  is  the 
stream  function  for  the  most  unstable  mode,  ^  is  that  of  the  first  subharmonic 
of  the  most  unstable  mode,  and  and  e ^  are  the  amplitudes  of  the  imposed  per¬ 
turbations.  The  presence  of  the  fundamental  mode  in  the  mixing  layer  produces 
a  single  vortex  rollup,  and  when  the  subharmonic  perturbation  is  superimposed, 
a  second  rollup  in  the  form  of  the  merging  of  two  vortices  occurs  (Ho  and 
Huerre,  1984).  The  vortex  rollup  and  pairing  have  dominant  effects  on  the 
chemical  reactions  and  the  rate  of  formation  of  the  products. 

The  downstream  boundary  is  an  artificial  computational  boundary  and  not  a 
natural  flow  boundary.  The  boundary  conditions  there  are  difficult  to  formu¬ 
late  in  a  rigorous  sense.  We  have  assumed  a  zero-gradient  condition  on  the 
velocity  and  species  fields.  This  is  not  an  exact  representation  of  the 
outflow  behavior.  However,  if  the  flow  velocity  is  outgoing  at  the  downstream 
boundary,  which  is  the  case  here,  the  vortlcity  and  species  transport  equations 
are  hyperbolic  when  the  molecular  diffusion  is  small.  Therefore,  the  upstream 
Influence  may  be  weak.  The  effects  of  this  outflow  boundary  condition  on  the 
solution  are  later  assessed  by  comparing  results  of  the  simulations  for  the 
same  flow  with  various  lengths  of  the  computational  domain. 

The  periodicity  of  the  flow  in  the  Y  direction  enables  us  to  employ  a 
pseudospectral  numerical  method  (Gottlieb  and  Orszag,  1977)  in  that  direction. 


In  the  s t Teamwise  direction,  however,  an  overall  second-order  finite  difference 
scheme  (Leonard,  1970)  is  employed  with  applications  of  a  quadratic  upwind 
differencing  for  the  convection  term  where  the  leading  dispersive  truncation 
error  of  the  second-order  central  scheme  is  removed.  This  scheme  is  expected 
to  be  accurate  for  the  fine  grids  employed  in  these  calculations.  We  use  the 
Adams-Bashforth  (Roache,  1972)  technique  for  time  discretization,  which  is  a 
second-order-accurate  scheme . 

The  computational  domain  was  selected  to  be  a  region  bounded  by  (0  <  X  <  32X) 
and  (-8X  <  Y  <  8X) .  There  are  128  equally  spaced  finite  difference  grid  points 
in  the  X  direction  and  64  Fourier  modes  in  the  Y  direction.  Since  the  details 
of  the  chemistry  are  neglected  and  a  fast  chemistry  model  has  been  employed, 
this  computational  domain  seems  to  be  accurate  enough  for  our  solution  proce¬ 
dure.  The  Incremental  time  step  was  selected  to  be  small  enough  to  ensure 
both  accuracy  and  stability.  The  value  of  At*  (At*  ■  At  AU/X)  was  selected  to 
be  0.0125. 

The  flow  is  conveniently  characterized  by  two  nondlmenslonal  parameters: 
first,  the  Reynolds  number.  Re  **  AU  v/X,  based  on  the  initial  Bhear  layer  thick¬ 
ness,  the  mean  velocity  difference  across  the  layer,  and  the  kinematic  visco¬ 
sity;  and  second,  the  velocity  ratio,  R  “  AU/CU^  +  U2).  The  Shvab-Zeldovich 
variable  is  characterized  by  the  value  of  the  Peclet  number,  Pe  ■  AU  D/X.  For 
the  calculations  reported  below,  we  have  selected  moderate  values  of  the 
Reynolds  and  Peclet  numbers  to  ensure  accuracy  of  the  numerical  simulations. 

As  discussed  by  Ho  and  Huerre  (1984),  the  two-dimensional  vortex  dynamics  of 
the  flow  is  not  very  sensitive  to  the  Reynolds  number,  and  only  the  shape  of 
the  vortices  depends  on  the  value  of  the  Reynolds  number  employed.  Therefore, 
simulations  even  at  moderate  Reynolds  numbers  provide  useful  information  about 
the  effects  of  the  dynamics  of  vorticity  on  chemical  reactions. 

3.  PRESENTATION  OF  RESULTS 
Effects  of  Coherent  Structures 

As  discussed  in  the  previous  section,  a  small  perturbation  at  the  inflow 
is  required  to  trigger  the  initial  vortex  rollup.  The  perturbation  levels, 
and  c^,  are  both  set  equal  to  0.05.  The  velocity  shear  parameter  R  is  equal 
to  0.5,  and  the  values  of  the  Reynolds  and  Peclet  numbers  are  set  equal  to  50. 


We  present  normalized  vorticity  contour  plots  at  t*  "  31  for  the  case  in 
which  the  most  unstable  mode  is  the  only  perturbation  applied  at  X  ■  0.  Figure  2 
indicates  the  formation  of  rolled-up  vortices  at  equal  wavelengths.  These 
vortices  diffuse  along  the  mixing  layer  as  indicated  by  the  decay  of  the  peak 
value  of  the  vorticity  downstream.  The  formation  of  the  rolled— up  vortices  is 
also  evident  in  the  contour  plot  of  the  Shvab-Zeldovich  variable  presented  in 
Figure  3.  From  these  two  figures,  it  is  clear  that  the  perturbation  at  the 
inlet  causes  the  rollup  of  the  vorticity  near  the  inlet.  As  the  fluid  is 
convected  downstream,  additional  vortices  are  created  at  equal  wavelengths 
from  each  other.  The  same  behavior  was  also  observed  in  numerical  and  experi¬ 
mental  studies  of  Mclnville  et  al.  (1985)  in  their  studies  of  the  large  co¬ 
herent  structures  in  a  forced  shear  layer.  As  the  vortices  reach  the  outflow 
boundary,  the  zero— gradient  condition  seems  to  allow  them  to  travel  out  of  the 
computational  domain.  Employment  of  this  weak  condition  at  the  outflow  may 
cause  some  errors  near  the  outflow  boundary.  Previous  numerical  experiments 
by  Givi  (1981),  however,  show  that  the  regions  of  the  error  associated  with 
this  condition  are  concentrated  at  a  small  area  near  the  outflow  and  do  not 
substantially  affect  the  inner  region  of  the  computational  domain.  To  assess 
the  validity  of  this  assumption,  some  computations  were  performed  in  two 
computational  domains.  The  transverse  dimensions  of  these  two  are  identical, 
while  one  of  the  cases  has  a  streamwise  domain  half  the  size  of  the  other. 

The  grid  spacing  in  the  streamwise  direction  is  maintained.  Comparison  of  the 
results  clearly  showed  that  the  boundary  conditions  imposed  at  the  outflow 
boundary  do  not  affect  the  solution  in  the  interior  of  the  computational 
domain.  Therefore,  it  seems  justified  that  a  zero-gradient  boundary  condition 
at  the  outflow  is  adequate  for  our  purposes. 

The  effects  of  adding  the  perturbation  corresponding  to  the  first  subhar¬ 
monic  mode  at  X  -  0  are  shown  in  Figure  4.  In  this  figure,  we  present  the 
rollup  and  merging  of  the  vortices  at  times  of  t*  ■  31  and  t*  ”  39.  The 
merging  of  the  vortices  caused  by  the  subharmonic  perturbation,  as  computed 
here,  has  been  observed  in  experiments  of  Ho  and  Huerre  (1984).  These  merged 
vortices,  again,  diffuse  and  smear  as  they  move  downstream,  and  no  further 
merging  is  observed.  The  same  behavior  was  also  observed  in  the  calculations 
of  Davis  and  Moore  (1985).  The  vortex  merging  associated  with  this  subharmonic 
mode  is  also  clear  from  the  Shvab-Zeldovich  variable  contour  plots  shown  in 


Figure  5.  Again,  the  errors  associated  with  the  outflow  boundary  are  concen¬ 
trated  only  near  the  outflow  and  do  not  seem  to  affect  the  interior  of  the 
computational  domain. 

Contour  plots  of  the  computed  product  concentrations  of  the  chemical 
reaction  are  presented  in  Figures  6  and  7.  In  Figure  6,  only  the  perturbation 
associated  with  the  most  unstable  mode  is  added  to  the  mean  flow,  and  in 
Figure  7,  the  perturbations  associated  with  the  fundamental  frequency  and  its 
first  subharmonic  are  added  together  to  the  mean  flow.  The  regions  of  high 
product  concentrations  are,  as  expected,  near  the  flame  sheet,  where  the  value 
of  the  Shvab-Zeldovich  variable  is  equal  to  its  stoichiometric  value  (f  “  0.5). 
The  vortex  rollup  brings  unreacted  species  together  from  both  streams  to  the 
chemical  reaction  zone.  This  region  is  marked  by  a  very  steep  gradient  of  the 
chemical  product  near  the  braids  of  the  vortices  but  is  more  diffusive  in  the 
core  of  the  vortex.  The  vortex  rollup  increases  the  interface  between  the  two 
streams  by  stretching  the  reaction  zone  and,  therefore,  increasing  the  amount 
of  products  generated.  A  comparison  of  these  two  figures  indicates  that  the 
vortex  merging  results  in  increased  formation  of  product  in  the  mixing  layer. 

In  these  computations,  the  increase  of  flame  sheet  surface  due  to  rollup 
of  vortices  is  about  97  percent  in  comparison  with  unforced  and  nondiffusive 
computations,  resulting  in  a  48-percent  increase  in  total  integrated  products 
formed  in  the  mixing  layer.  Merging  of  the  vortices  Increases  the  flame  sheet 
surface  area  by  163  percent  and  the  total  amount  of  products  by  69  percent  in 
comparison  with  the  unforced  and  nondiffusive  case.  Bearing  in  mind  the  low 
level  of  initial  disturbance,  this  high  level  of  response  may  one  day  be 
exploited  in  a  practical  device. 

Comparison  of  Figures  6  and  2  and  also  Figures  7  and  4  shows  the  similarity 
between  the  profiles  of  the  vorticity  and  the  chemical  product.  Both  are  high 
near  the  reaction  zone,  gradually  diffusing  outward  with  very  similar  profiles. 
This  is  to  be  expected,  since  the  chemical  reaction  is  to  occur  where  the  two 
reactants  meet  near  the  interface,  which  is  also  a  region  of  high  velocity  gra¬ 
dient.  The  same  conclusions  were  also  drawn  from  earlier  work  on  spectral  cal¬ 
culations  of  temporally  growing  mixing  layers  (Riley  et  al.,  1986).  However, 
the  asymmetric  nature  of  the  mixing  mechanism  in  the  shear  layer  is  well  repre¬ 
sented  in  the  spatially  developing  calculations,  which  are  discussed  next. 


Statistical  Variables 

The  results  of  direct  numerical  simulations  have  been  studied  to  examine 
the  validity  of  some  of  the  previously  used  turbulence  models  in  which  the 
large  coherent  structures  have  not  been  taken  into  account.  In  previous  work 
(Givi  et  al. ,  1984,  1985),  many  different  possibilities  for  the  closure  of  the 
hydrochemical  equations  in  parallel  shear  flows  were  discussed.  Among  the 
models  considered,  the  one  with  the  best  overall  performance  was  the  combina¬ 
tion  of  a  modeled  joint  probability  density  function  (pdf)  for  the  scalar 
variables  and  the  K-e  model  of  turbulence  (Launder  and  Spalding,  1972)  for 
the  hydrodynamic  closure.  A  formulation  based  on  a  transport  equation  for  the 
pdf  of  scalar  variables  has  the  advantage  that  the  effects  of  chemical  reac¬ 
tions  appear  in  a  closed  form.  However,  models  are  needed  for  the  closure  of 
the  molecular  mixing  term  and  also  the  turbulent  convection.  Givi  et  al. 

(1985)  employed  a  qualesence/dispersion  model  for  the  closure  of  the  molecular 
mixing  term  and  a  gradient  diffusion  approximation  for  the  closure  of  the 
turbulent  flux  of  the  pdf.  The  results  were  compared  with  experimental  data 
of  Masutanl  and  Bowman  (1984)  in  similar  hydrochemical  conditions. 

As  far  as  the  first  few  moments  are  concerned,  the  results  obtained  from 
the  pdf  transport  equation  are  in  reasonable  agreement  with  experimental  data. 
However,  a  major  difference  between  the  calculated  and  measured  Bhapes  of  the 
probability  density  function  is  observed.  This  is  indicated  in  Figure  8.  In 
this  figure,  the  pdf’s  of  a  nonreacting  scalar  variable  (similar  to  the  varia¬ 
ble  f  defined  above),  are  presented  at  a  streamwlse  location.  The  profile 
obtained  from  experimental  measurements  of  Masutanl  and  Bowman  (1984)  is  also 
presented  on  the  same  figure.  The  experimental  results  Indicate  an  interme¬ 
diate  peak  at  a  fixed  value  of  the  normalized  concentration  between  the  two 
delta  functions  at  0  and  1.  In  the  predicted  pdf's  two  spikes  also  exist  at 
the  normalized  concentrations  of  0  and  1,  Indicating  the  concentrations  at 
freestreams.  The  height  of  the  pdf  in  the  middle  region  is  in  reasonably 
correct  order-of-magnitude  agreement  with  the  experimental  data.  However,  the 
location  of  the  predicted  pdf  peak,  with  respect  to  the  normalized  concentra¬ 
tion,  gradually  shifts  as  the  layer  is  traversed.  This  is  a  major  difference 
between  the  predicted  and  measured  pdf's,  as  the  modeled  pdf  transport  equation 
is  unable  to  predict  the  bimodal  pdf's  observed  experimentally.  As  suggested 
by  Givi  et  al.  (1985),  Koochesfahani  (1983),  and  Masutanl  (1985),  this  discre¬ 
pancy  is  mainly  due  to  the  shortcomings  associated  with  the  gradient  diffusion 


modeling  of  the  turbulent  flux  of  the  pdf.  In  highly  intermittent  flows  such 
as  mixing  layers,  the  continuity  of  turbulent  flow  is  interrupted  by  the 
presence  of  the  nonturbulent  surrounding  flows.  Therefore,  a  simple  gradient 
diffusion  model  is  not  expected  to  account  accurately  for  this  discontinuity. 

Direct  comparison  of  the  shape  of  the  pdf  calculated  from  direct  numerical 
simulations  with  that  of  experimental  measurements  is  not  possible  due  to  the 
different  Reynolds  numbers  used.  However,  a  great  deal  of  understanding  of  the 
Influence  of  the  large-scale  structures  in  the  mixing  mechanism  of  the  shear 
layer  can  be  gained  by  examining  the  computed  concentration  pdf.  The  profile 
of  the  pdf's  calculated  from  the  results  of  the  direct  numerical  simulation  at 
the  streamwise  location  of  X  ■  24X  is  presented  in  Figure  9.  This  location 
corresponds  to  the  point  where  the  first  merging  between  the  two  neighboring 
vortices  occurs.  It  is  Indicated  in  this  figure  that  at  any  cross-stream 
location  away  fr^m  the  freestreams,  the  pdf  has  approximately  three  peak 
values.  The  first  peak  is  closer  to  the  low-speed  fluid  concentration,  the 
second  peak  is  closer  to  the  high-speed  fluid  concentration,  and  the  third 
peak  is  in  a  mixed  concentration  between  the  other  two  but  closer  to  the 
high-speed  side.  The  value  of  the  pdf  of  the  mixed  concentration  is  larger 
than  the  pdf  corresponding  to  other  concentrations;  also,  the  pdf  of  the 
concentration  closer  to  the  high-speed  fluid  side  has  a  larger  peak  than  the 
pdf  of  the  low-speed  concentration.  This  trimodal  nature  of  the  pdf,  and  also 
the  fact  that  the  mixed  fluid  pdf  has  a  concentration  near  the  high-Bpeed  side 
concentration,  indicates  the  asymmetry  of  entrainment  due  to  large-scale 
structures  in  the  mixing  region  of  the  shear  layer  and  also  Indicates  that 
more  fluid  from  the  high-speed  stream  than  the  low-speed  side  is  drawn  into 
the  mixing  core  of  the  layer.  Comparison  of  this  figure  with  Figure  8  indi¬ 
cates  that  the  pdf's  obtained  here  have  a  structure  closer  to  the  experimental 
observations  of  Hasutani  and  Bowman  (1984)  than  the  previous  calculations 
employing  a  modeled  joint  pdf  transport  equation  (Givi  et  al.,  1985).  However, 
the  calculated  predominant  peak  of  the  pdf  at  the  mixed  fluid  concentration  is 
less  pronounced  than  that  observed  experimentally  by  Masutani  and  Bowman 
(1984).  This  is  possibly  due  to  the  fact  that,  in  the  present  two-dimensional 
calculations,  only  the  isolated  effects  of  the  large  coherent  structures  are 
examined,  and  the  effects  of  three-dimensional  random  turbulent  fluctuations 
are  not  presented.  Nevertheless,  the  results  obtained  from  direct  numerical 
simulations  indicate  that  the  mixing  mechanism  in  the  core  of  the  shear  layer 


is  asymmetric,  and  the  mixed  fluid  has  a  concentration  closer  to  the  high-speed 
fluid  concentration.  This  agrees  with  the  experimental  observations  of 
Masutani  and  Bowman  (1984)  and  Koochesfahani  (1983).  Three-dimensional  flow 
simulations  with  random  turbulent  fluctuations  are  currently  under  study  and 
are  expected  to  result  in  better  agreement  with  experimental  data. 

The  mean  and  rms  values  of  the  conserved  scalar  obtained  from  the  inte¬ 
gration  of  the  pdf  profiles,  presented  above,  are  shown  in  Figure  10.  This 
figure  shows  that  the  mean  value  of  the  concentration  has  an  apparent  plateau 
in  the  middle  region  of  the  shear  layer.  This  "plateau"  is  at  a  concentration 
where  a  predominant  peak  occurs  in  the  pdf  (Figure  9).  Similar  behavior  has 
also  been  observed  in  experimental  measurements  (Koochesfahani,  1983)  and 
indicates  the  influence  of  the  large-scale  structure  on  the  mean  value  of  the 
concentration.  The  resulting  rms  values  of  the  normalized  non— reacting  con¬ 
centrations  across  the  layer,  obtained  from  the  direct  numerical  simulations, 
indicate  that  there  are  two  apparent  maximums  in  the  rms  profile.  These  two 
points  correspond  to  the  location  where  the  gradient  of  the  mean  value  is 
highest.  The  same  behavior  was  also  observed  in  the  experimental  results  of 
Masutani  (1985).  Calculations  using  a  gradient  diffusion  model  for  the  pdf 
turbulent  transport  are  unable  to  predict  this  detailed  behavior,  as  the 
calculations  of  Givi  et  al.  (1985)  indicate.  The  predicted  results  of  the 
mean  concentration  by  Givi  et  al.  (1985)  form  a  fairly  smooth  bell-shaped 
profile  with  a  much  less  clear  double  "hump"  in  the  middle  region.  However, 
the  pdf  calculation  does  show  that  the  location  of  the  rms  maximums  coincides 
with  the  location  of  the  highest  mean  value. 

From  the  discussion  given  above,  it  is  clear  that  the  intermittency  caused 
by  the  large  coherent  structures  contributes  greatly  to  various  vital  statis¬ 
tics  of  turbulent  flows.  The  results  of  direct  numerical  simulations  can 
illuminate  the  qualitative  behavior  of  large-scale  structures  in  intermittent 
flow,  such  as  the  spatially  developing  mixing  layer  considered  here.  Quantita¬ 
tive  comparison  with  the  experimental  data,  however,  is  not  yet  possible.  A 
three-dimensional  flow  simulation  may  be  required  for  meaningful  quantitative 
comparison  with  experimental  results. 
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4.  CONCLUSIONS  AND  FUTURE  WORK 

Numerical  simulations  have  been  performed  for  the  large-scale  vortex 
dynamics  in  a  forced  spatially  developing,  reacting  mixing  layer.  A  combina¬ 
tion  of  a  pseudospectral  and  a  second-order  finite  difference  technique  has 
been  employed  to  study  the  effects  of  the  rollup  of  the  vortices  and  of  their 
merging  on  the  rate  of  formation  of  the  reaction  product.  As  an  earlier 
simulation  (Riley  et  al. ,  1986;  McMurtry  et  al.,  1985)  using  a  temporally 
developing  approximation  indicates,  the  vortex  rollup  and  merging  cause  the 
stretching  of  the  interface  between  the  reactants,  and  thereby  increase  the 
stream  surface  area  on  which  the  chemical  reaction  occurs.  In  addition,  the 
present  simulations  of  a  realistic  spatially  developing,  chemically  reacting 
mixing  layer  exhibit  the  asymmetric  properties  of  the  entrainment,  which  have 
been  observed  in  experimental  results  (Masutani  and  Bowman,  1984). 

The  results  of  direct  numerical  simulations  reported  here  and  recent 
experimental  measurements  Indicate  the  importance  of  large  coherent  structures 
in  the  mixing  region  of  the  shear  layer.  This  suggests  the  need  for  better 
turbulence  models  in  order  to  predict  accurately  the  mechanism  of  mixing  and 
chemical  reaction  in  intermittent  flows  such  as  mixing  layers. 

Much  remains  to  be  done  in  our  continuing  effort  in  numerical  combustion 
simulations.  Consideration  of  a  temperature-dependent  reaction  rate  with 
intermediate  Damkohler  numbers  has  just  been  successfully  finished  for  a 
temporally  growing  mixing  layer  (Givi  et  al. ,  1986)  and  should  be  included 
with  the  spatially  growing  mixing  layer  studied  in  thiB  paper.  This  simula¬ 
tion  can  address  the  questions  regarding  the  quenching  and  extinction  of 
diffusion  flames.  Next,  three-dimensional  calculations  Involving  a  shear 
layer  with  a  random  turbulence  field  will  be  studied.  This  is  a  step  toward 
simulations  of  turbulent  flow  and  will  provide  a  better  understanding  of  the 
physics  of  turbulent  diffusion  flames.  The  effects  of  chemical  heat  release, 
at  a  temperature-dependent  reaction  rate,  on  the  flow  field  and  product  forma¬ 
tion  in  a  spatially  developing  mixing  layer  would  also  be  very  interesting  and 
will  be  addressed  in  future  papers.  Finally,  it  would  be  interesting  to 
examine  the  effects  of  the  vortex  merging  further  by  looking  at  the  influence 
of  even  higher  subharmonics  and  also  the  effects  of  phase  relations  between 
the  different  modes  on  the  chemical  reactions  and  the  flame  convolution  in  a 
spatially  evolving  mixing  layer. 
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